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Abstract. Centered numerical fluxes can be constructed for compressible Euler equa- 
tions which preserve kinetic energy in the semi-discrete finite volume scheme. The es- 
sential feature is that the momentum flux should be of the form f m , = p. , \ +u , ■_, 1 f p , 

where u- + \ — («y+wy +1 )/2 and pj + i,f P - + \ are any consistent approximations to the 

pressure and the mass flux. This scheme thus leaves most terms in the numerical 
flux unspecified and various authors have used simple averaging. Here we enforce 
approximate or exact entropy consistency which leads to a unique choice of all the 
terms in the numerical fluxes. As a consequence novel entropy conservative flux that 
also preserves kinetic energy for the semi-discrete finite volume scheme has been pro- 
posed. These fluxes are centered and some dissipation has to be added if shocks are 
present or if the mesh is coarse. We construct scalar artificial dissipation terms which 
are kinetic energy stable and satisfy approximate /exact entropy condition. Secondly 
we use entropy-variable based matrix dissipation flux which leads to kinetic energy 
and entropy stable schemes. These schemes are shown to be free of entropy violating 
solutions unlike the original Roe scheme. For hypersonic flows a blended scheme is 
proposed which gives carbuncle free solutions for blunt body flows. Numerical results 
for Euler and Navier-Stokes equations are presented to demonstrate the performance 
of the different schemes. 
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1 Introduction 

The numerical solution of compressible Euler and Navier-Stokes (NS) equations by the 
finite volume method is now a routine task in many industries. Due to their non-linear 
hyperbolic nature, solutions of Euler equations can be discontinuous with the presence 
of shocks or contact discontinuities. Discontinuous solutions must necessarily satisfy the 
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Rankine-Hugoniot jump conditions which are a consequence of conservation. However 
it is well known that discontinuous solutions that satisfy the jump conditions can be still 
non-unique and an additional entropy condition has to be imposed in order to select the 
unique weak solution. In the case of Euler equations, there is a natural entropy condition 
which comes from the entropy condition in thermodynamics which must also be satisfied 
by the numerical scheme. Additionally other global balance equations like that for the 
total kinetic energy must also be consistently approximated by the numerical solutions. 
The finite volume method requires the computation of the inviscid and viscous fluxes 
across the boundaries of the finite volumes. The design of these fluxes must incorpo- 
rate the properties of the Euler/NS equations like entropy condition and kinetic energy 
preservation. There exists a vast library of numerical flux functions for the Euler equa- 
tions and some of these like the Godunov scheme and kinetic scheme can be shown to 
satisfy the entropy condition. The popular Roe scheme |20] does not satisfy the entropy 
condition and can give rise to entropy violating shocks near sonic points. Various en- 
tropy fixes for Roe scheme have been proposed which involve preventing the numerical 
dissipation from vanishing at sonic points. Tadmor |25[ proposed the idea of entropy 
conservative numerical fluxes which can then be combined with some dissipation terms 
using entropy variables to obtain a scheme that respects the entropy condition, i.e., the 
scheme must produce entropy in accordance with the second law of thermodynamics. 
This is a more mathematically rigorous approach to construct entropy stable schemes 
for conservation laws. However some of these entropy conservative numerical fluxes 
have to be computed with quadrature rules since the integrals involved in the definition 
of the flux cannot be evaluated explicitly. For the Euler equations, Roe proposed explicit 
entropy conservative numerical fluxes | lOpl) which are augmented by Roe-type dissipa- 
tion terms using entropy variables. These schemes do not suffer from entropy violating 
solutions that are observed in the original Roe scheme. However for strong shocks, even 
the first order schemes can produce oscillations indicating that the amount of numerical 
dissipation is not sufficient. Roe 1 10 1 proposed modifying the eigenvalues of the dissipa- 
tion matrix which lead to non-oscillatory solutions. The modification of the eigenvalues 
is such that the amount of entropy production is of the correct order of magnitude for 
weak shocks. The availability of cheap entropy conservative fluxes allows us to use the 
procedure of |15| to develop high order accurate entropy conservative schemes. Matrix 
dissipation can be added following the ENO procedure of |4j to develop arbitrarily high 
order accurate entropy stable schemes for the Euler equations on structured grids. 

Faithful representation of kinetic energy evolution is another desirable property of a 
numerical scheme |l2| . This is important for direct numerical simulation of turbulent 
flows where the kinetic energy balance plays an important role in the evolution of tur- 
bulence |16}18}22| . The scheme is also stable in the sense that spurious kinetic energy is 
not produced by the numerical fluxes. Godunov schemes are found to have wrong order 
of kinetic energy dissipation and entropy production which leads to excessive damping 
of flow structures p7| . One way to construct kinetic energy preserving schemes is to use 
a skew symmetric form for the non-linear convective terms combined with finite differ- 
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ence operators which satisfy summation by parts property [3,8,14]. The skew symmetric 
form is also found to yield smaller aliasing errors which makes them more robust [3|. 
For compressible simulations, it is necessary to also satisfy conservation of mass, mo- 
mentum and energy and hence it is preferable to solve the Navier-Stokes equations in 
conservation form using the finite volume method. The essential feature for a numeri- 
cal flux in a semi-discrete finite volume method to correctly capture the kinetic energy 
balance is that the momentum flux should be of the form f m , = p- . 1 +w , , 1 f p , where 

Uj + l = (uj + Uj + i)/2 and Pj + i,f P + i are any consistent approximations to the pressure and 

the mass flux. This scheme thus leaves most terms in the numerical flux unspecified and 
various authors have used simple averaging. Subbareddy and Candler p4| have pro- 
posed a fully discrete finite volume scheme for the compressible Euler equations which 
preserves kinetic energy but the resulting scheme is implicit. All of these kinetic energy 
preserving schemes are however not entropy conservative, while on the other hand, the 
entropy conservative schemes do not have the kinetic energy preservation property. It is 
thought that for DNS of compressible flows, a numerical scheme which preserves kinetic 
energy and satisfies entropy condition is desirable since such schemes would be non- 
linearly stable. Schemes which satisfy entropy condition are found to lead to stable den- 
sity fluctuations in compressible isotropic turbulence simulations, while schemes which 
do not have this property can be unstable with respect to density fluctuations )8p8"|. Us- 
ing a canonical splitting of the flux with exponential entropy, Gerritsen and Olsson [5| 
construct entropy stable schemes which are however not conservative with respect to 
mass, momentum and energy. In J8j, skew-symmetric form of convective terms is used 
to enforce better entropy consistency which is then shown to lead to schemes capable of 
computing high Reynolds number turbulence without artificial dissipation or filtering. 

In the present work, we construct explicit centered numerical fluxes for the compress- 
ible Euler equations which are entropy conservative and also preserve kinetic energy in 
the case of the semi-discrete finite volume scheme. All the numerical fluxes presented 
here can be used in any finite volume scheme on structured or unstructured grids. Two 
versions of the numerical flux are constructed, one of which is only approximately en- 
tropy consistent but has simpler expressions, while the second one is exactly entropy 
conservative but involves certain logarithmic averages requiring more computations. 
Due to lack of upwinding, the schemes are not stable for discontinuous solutions and 
for NS equations on coarse meshes for which shocks may not be well resolved. They 
yield stable solutions for Navier-Stokes equations when used on very fine meshes where 
the physical viscosity is enough to stabilize the scheme. However for Euler equations 
and for NS equations on coarse meshes, the centered fluxes are unstable and must be 
augmented with dissipation terms. Firstly, we construct scalar dissipation terms using 
second and fourth order differences as in the JST scheme fll3) . The second order dissipa- 
tion terms which are active near shocks are kinetic energy and entropy stable while the 
fourth order dissipation is active only in smooth regions of the flow. Secondly, we also 
use entropy variable based matrix dissipation flux similar to the Roe scheme which can 
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be shown to lead to entropy generation [25|. The eigenvalue modification of Roe |T0| is 
used to compute strong shocks without oscillations. All the schemes are shown to give 
entropy consistent solutions in cases where the Roe scheme would give entropy violating 
shocks. The matrix dissipation scheme can also be made kinetic energy dissipative by a 
proper choice of the eigenvalues appearing in the dissipation matrix. The entropy stable 
schemes with matrix dissipation preserve stationary contacts exactly but also suffer from 
1-D shock instability and the carbuncle phenomenon. A modification of the eigenvalues 
in the dissipation flux is suggested which avoids these problems but is still able to ac- 
curately compute shear flows like boundary layers. The performance of the schemes is 
demonstrated on inviscid and viscous test cases involving shocks and contact disconti- 
nuities. 

The rest of the paper is organized as follows. Section (p) introduces the 1-D Navier- 
Stokes equations and finite volume method and section (pT discusses the kinetic energy 
preservation property. The entropy condition is introduced and the entropy conserva- 
tive fluxes are derived in section Q together with some justification for their unique- 
ness. Scalar artificial dissipation terms which satisfy kinetic energy and entropy sta- 
bility are derived in section ^ while matrix dissipation flux is treated in section ([6]!. 
Section j7]) shows numerical results for stationary 1-D shock problem and modifications 
of the scheme for monotone resolution of shocks, while section ^ introduces a hybrid 
scheme. The 2-D equations and numerical fluxes are discussed in section and sec- 
tion pD) ends with some numerical results on a range of test problems. 



2 1-D NS equations and finite volume scheme 

The one dimensional Navier-Stokes equations can be written in vector conservation form 
as 

dt dx dx 



(2.1) 



where u is the set of conserved variables and f, g are the inviscid and viscous fluxes 
whose expressions are given by 
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In the above equations p is the density, u is the velocity, p is the pressure and E is the 
total energy per unit volume which for a perfect gas is given by E = p/ { r y — \)+pu 2 /2, 
where 7 is the ratio of specific heats which is taken to be constant. Moreover, t, q are the 
shear stress and heat flux for which we take the Newtonian and Fourier laws respectively, 
leading to 

4 du dT 
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where ]i, k are the coefficient dynamic viscosity and heat conduction respectively. In the 
absence of the viscous fluxes g, the resulting Euler equations form a hyperbolic system 
of conservation laws. Consider a partition of the one dimensional domain into uniform 
finite volumes of size Ax and the /'th cell is the interval (x- i.x - , 1). The semi-discrete 

7 2 • ' 2 

finite volume scheme is given by 

Ax ^ +f /+r f H=s, + rsH (Z3) 

where u, is the cell average value in the /'th cell and f ■ , i. e, , i are numerical inviscid 

J o ■< ]-r 2 °/+2 

and viscous fluxes respectively at the interface In the numerical computations, 

the above set of ordinary differential equations will be solved using a strong stability 
preserving Runge-Kutta scheme [23| . 



3 Kinetic energy preserving scheme 



The kinetic energy is an important quantity in fluid flows and it is destroyed by the phys- 
ical viscosity. In turbulent flows, the kinetic energy injected into the fluid at large scales 
cascades to smaller scales and is eventually destroyed by viscosity. Hence it is desirable 
that the numerical scheme faithfully represent the kinetic energy balance consistent with 
the Navier-Stokes equations. The kinetic energy per unit volume K= \pu 2 satisfies the 
following equation 
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where we have ignored boundary conditions. The first term on the right is the rate at 
which work is done by pressure forces and is present only for compressible flows. The 
second term represents the irreversible destruction of kinetic energy which is converted 
into internal energy due to viscous dissipation. We would like the numerical scheme to 
also satisfy this equation in a discrete sense which will then be refered to as kinetic energy 
preserving scheme. Consider the following approximation for the inviscid and viscous 
fluxes 
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Throughout the paper, we will use the overbar to denote the arithmetic average. The 
quantities p,fP,u,f e are assumed to be consistent approximations but are yet to be speci- 
fied. The kinetic energy varies in time according to 
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We now derive the global kinetic energy balance equation for the finite volume scheme 
by adding the equation from each cell and summing over all cells in the grid. 




This is consistent with the continuous kinetic energy equation. The crucial property used 
in the above proof was that the momentum flux has the form f m = p+ufP which leads to 
the disappearance of the convective flux from the kinetic energy equation. However, we 
still have freedom in the choice of fP,p,u,f e . Jameson [12| makes the following choice: 

f p = pu, p = p, f = puH 

where H is the enthalpy. We will refer to this as the KEP flux. In the present work we will 
determine the numerical flux in a unique manner from entropy considerations. 



4 Entropy condition 

The concept of entropy condition is borrowed from the second law of thermodynamics 
and generalized to any arbitrary system of hyperbolic conservation laws [7[. Let U(u) be 
a strictly convex function; then U(u), F(u) is said to be an entropy-entropy flux pair if 
they satisfy 

U'(u)f'(u) = f'(u) 

Then smooth solutions of the inviscid equations satisfy an additional conservation law 

dll dF n 

But for discontinuous solutions we can only satisfy the entropy inequality 
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where the inequality is satisfied in a weak sense. Define the entropy variables as 

v(u) = U'(u) 

and since U(u) is strictly convex, we can invert the above equation to obtain u = u(v). 
Define the quantity i/>(v) which is the dual of the entropy flux F(u) by the relation 

<Kv)=vf(u(v))-F(u(v)) 

The finite volume scheme for the hyperbolic conservation law is given by 

du,- 

AX— i +f;, 1 -f, 1=0 

Taking the dot product of the above equation with v,-, we obtain the entropy equation 

Tadmor \ 25 26) introduced the idea of an entropy conservative numerical flux which should 
satisfy the following condition 



Then we can write 



where 
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is a consistent numerical entropy flux. The semi-discrete scheme satisfies the entropy 
conservation equation 

dUj 

Ax— } -+F j+1 -F i _i=0 
at i + 2 j 2 

An entropy conservative flux is given by Tadmor p5| as 
f 



M = l! i( ~ V j^ d6 ' v /+| (0)=v / +0(v /+1 -v / ) 



which usually requires some numerical quadrature to approximate the flux. Once an 
entropy conservative flux has been constructed, we can add dissipative terms to the flux 
which leads to the satisfaction of the entropy condition as given by equation (4.1 1, i.e., the 
dissipative flux must lead to generation of entropy. 
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4.1 Entropy condition for Euler equations 

For the Euler equations we can take the entropy-entropy flux pair to be 

u=--eV f=-^t (4.3) 

7 — 1 7 — 1 

where s is the physical entropy given by 

s = ln(p)-ln( (D T)+const = -(7-l)ln(p)-ln( j 6)+const / P = ^f ( 44 ) 

and the constant term can be ignored. There are many other possible choices for the 
entropy function U(u) = prj{s) where tj is any convex function |7j, but the above choice 
is the only one which is consistent with the entropy condition from thermodynamics (9j 
in the presence of heat transfer. Since we work with the correct choice of the entropy 
function, the schemes we develop will satisfy entropy condition for the Navier-Stokes 
equations also. The entropy variables v and the Legendre transform xp are given by 



2f>u 



ip = m = pu (4.5) 



Hence an entropy conservative numerical flux for the Euler equations has to satisfy the 
following condition 

(v 7 +i-v ; -)-f - + i =m j+1 -mj (4.6) 

This provides only one equation whereas the flux f has more than one component; we 
can expect that there are many possible entropy conservative fluxes. 



4.2 Roe's entropy conservative flux 



Roe 1 21 J has constructed explicit numerical fluxes for the Euler equations which satisfy 
condition (4.6 1 for the entropy given in equation ( |4.3| but do not have the kinetic energy 
preservation property. Introduce the set of independent state variables or parameter vec- 
tor 
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Roe also introduces the logarithmic average <p of any strictly positive quantity q>i, (p r 
which is defined as 

<p r -<pi = Af 
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A numerically stable procedure to compute the average when cpi cp r is given in [10| 
Then the entropy conservative numerical flux at any cell face ] + \ is given by 
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and all the averages are evaluated using the state in cell j and 7+I. This flux is entropy 
conservative and it can be made entropy stable [21 1 by adding entropy variable based 
matrix dissipation terms which is described in a later section. The resulting numeri- 
cal flux will be refered to as ROE-ES flux. Note that the momentum flux contains the 
weighted average velocity u while kinetic energy preservation requires the presence of 
the arithmetic average u. Hence the above entropy conservative flux is not kinetic energy 
preserving. In the following sections, we construct numerical fluxes by approximately or 
exactly satisfying the entropy condition for the Euler equations which also preserves the 
kinetic energy. 



4.3 Derivation of approximately entropy consistent flux 

We now derive kinetic energy preserving flux which is approximately entropy consistent. 
In order to simplify the notations, we will drop the subscripts on the cell indices, and use 
the convention that A ( • ) = ( • — ( • )y denotes the jump across the cell face /+ \. Equation 
<|4.6|l can then be written as 
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We can easily check that the following exact linearizations are true. 

A(j3u)=/3Aw+T/A^, A(pu)=pAu + uAp, A(/3m 2 ) =^Am 2 + m^Aj6 = 2j6mAu + i?A^ 



The entropy difference is 

As = -(7-l)Aln l o-Aln 1 6 
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Assuming smooth solutions we approximate the above differences by a Taylor formula, 
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and similarly for Aln/S. We then obtain 
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For smooth solutions, the higher order terms are of 0(Ax 3 ). Neglecting the third order 
terms, we try to satisfy condition (|4.6| which becomes 
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Since p,u,f> are independent variables, the above equation is satisfied if we choose 
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With the above unique choice, the numerical flux satisfies condition ( |4.6| to third order 
accuracy, i.e., 

(v /+1 -v ; )-f /+ i= m/+1 - m; +0(Ap) 3 +0(A/3) 3 

It is easy to check that the above numerical fluxes are consistent. We remark that these 
flux formulae have simpler expressions since they make use of arithmetic averages while 
the entropy conservative fluxes have slightly more complicated expressions with loga- 
rithmic averages as we see from Roe's fluxes and also the new flux derived in the next 
section. 
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4.4 Kinetic energy preserving and entropy conservative flux 

In this section we derive numerical fluxes which preserves the kinetic energy and exactly 
conserves the entropy. The jump in the entropy variables can be written as 



Av x = Aln(,0) + — j-Aln^-A^w 2 

Av 2 = 2 i 6Au+2uA i 6 
Av 3 = -2AjS 
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where we have made use of the logarithmic averages p and p as introduced in sec- 
tion (|4.2[). The flux is entropy conservative if equation (|4.6|l is satisfied, i.e., if 



fAvt +f m Av 2 +f e Av 3 = A(pu) = pAu + uAp 



(4.7) 



As in the previous section, the above equation is a finite difference equation involving 
differences in the independent variables p,u,f>. Equating the terms containing Ap on both 
sides yields the mass flux as 

fP = pu 

Next equating the terms containing Au we get the momentum flux as 
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The above form of the momentum flux satisfies the requirement for kinetic energy preser- 
vation. Finally, from the terms containing Aj6, we get the energy flux 
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These fluxes are consistent and have exactly the same form as the approximately con- 
sistent fluxes derived in the previous section except that the logarithmic averages p, f> 
are used instead of the arithmetic averages p, f> in the mass and energy fluxes. We can 
also see that the new entropy conservative flux is computationally less expensive as com- 
pared to Roe's entropy conservative flux given in section ( |4.2[ | since it does not require a 
parameter vector in its definition. 



4.5 Uniqueness of the numerical flux 

Roe 1 21 1 has proposed explicit entropy conservative flux for Euler equations and in the 
present work we have proposed another one. Tadmor |25) gives a recipe for the construc- 
tion of entropy conservative fluxes for any hyperbolic system endowed with an entropy- 
entropy flux pair. However, the new flux function proposed here also satisfies kinetic 
energy preservation. Since there are many numerical flux functions that satisfy entropy 
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conservation, it is interesting to ask if there is a unique numerical flux function that satis- 
fies both conservation properties and here we give some justification in favour of this for 
the case of two point numerical fluxes. The crucial property for kinetic energy preserva- 
tion is that the momentum flux must be of the form f' n = p+ufP where p is any consistent 
approximation to the pre ssure and fP is any consistent approximation to the mass flux. 
Looking at equation (4.71, we see that the term ufP can only come from the term A(fiu 2 ) 
in Avi. In the above derivation, we have written A(fiu 2 ) = u 2 Afi+fiAu 2 = u 2 Afi+2fiuAu 
which gives the desired form in the momentum flux but there are other ways to do the 
linearization. If instead we write this term as 



A(fiu 2 ) = fiuAu+uA(fiu) = {fiu + ufi)Au + u 2 Afi 

then we do not get the correct form of the momentum flux. Roe writes it in the form 
A(/3u 2 ) =A(y / J 8u) 2 = 2y / /3uA(y / /3u) which also does not have the correct form for kinetic 
energy preservation. 

The other type of non-uniqueness can arise with the choice of independent state vari- 
ables. In the derivation of the new flux, we wrote condition ( |4.6| in terms of jumps in p, 
u, fi. The other possibility is to use the variables p, u, p and we show that this does not 
lead to a proper definition of the numerical fluxes. The jump in the physical entropy at a 
cell face - 



2 can be written as 
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Then the jump in the entropy variables is given by 
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If we now use condition ( |4.6| , then the terms containing Au lead to the momentum flux 

fn = ]L + ufP 
J 2fi J 
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which satisfies the condition for kinetic energy preservation. The terms containing Ap 
yield the energy flux 

PjPj+i 1 : 



while the terms containing Ap yield 

1 



fP + uf n 
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From the last two equations and using the relation f = p/ (p/py+i), we obtain the mass 
flux as 

pI7 



(7-1) (y-l)pp 



The above mass flux depends on the ratio of specific heats 7 which is not physically 
meaningful. 

We can also use the variables p,u,fi as independent variables. In this case we have the 
following exact linearizations 

A01 = / ^ + ^— / ^-2TiMu 
V 7-1 jS 

Av 2 = 2(I/A0+j8Au) 

Aj7 3 = -ZA^ 

A(pu) = pAu+2u pAp+lufiAp 



Then condition (|4.6| gives the fluxes as 
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While this has the kinetic energy preservation property for the momentum flux, we find 
that the energy flux is not consistent. This indicates that the variables p, u, are the 
appropriate set for the satisfaction of the entropy conservative flux condition ( |4.6| for the 
Euler equations. 



4.6 Entropy equation 



Now we can derive the global entropy balance equation for the semi-discrete finite vol- 
ume scheme for the Navier-Stokes equations. Taking dot product of equation ( |2.3| with 
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the entropy variables Vj and summing up over all the finite volumes yields 



du 
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and the O(Ax) 3 terms are not present if the entropy conservative flux is used. The terms 
involving F- + i represent convection of entropy and cancel one another when we sum 
over all cells. The terms on the right which consist of viscous shear stress and heat flux 
can be shown to lead to entropy generation. 
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Hence the entropy equation becomes 
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If the entropy conservative flux is used, the O(Ax) 3 terms are not present and the entropy 
condition is satisfied exactly for any mesh size Ax, which is consistent with the entropy 
condition from the second law of thermodynamics. 

4.7 Summary of flux formulae 

We now list the approximately and exactly entropy consistent centered fluxes for the 1-D 
Euler equations. 



(1) The centered, kinetic energy preserving and entropy consistent numerical flux which 
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will be denoted by f * is given by 
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r*,0 , — r- 



*,m 
l 

2 



(4.8) 
(4.9) 

(4.10) 



where 



Pj+i 



This is equivalent to using the harmonic average for the temperature, i.e., 



(4.11) 



Pj + l= R Pj + \Tj 



]+ 



_ 27)7) +1 
- ^ T j + T j+1 



(2) The centered, kinetic energy preserving and entropy conservative numerical flux is 
given by 



r*,m 
r*,e 



. — r *,p 



f*'P _L77 f*' m 



(4.12) 
(4.13) 

(4.14) 



where pj + i is given by equation (4.11 1. 



5 Scalar artificial dissipation flux 

The entropy consistent fluxes have a central character and will lead to a stable scheme 
for Navier-Stokes equations only on highly resolved meshes. For the scheme to be stable 
on coarse meshes and inviscid problems, let us introduce second order, scalar artificial 
dissipation term into the inviscid fluxes, so that the numerical flux becomes 

Vi=^r d=[dp ' Dm - De]T - A -° 

We will derive the form of the dissipation from kinetic energy and entropy stability con- 
soderations. 
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5.1 Kinetic energy stability 

We make the following obvious choice for the dissipation in mass and momentum fluxes 

DP = Ap, D m = A(pu) 



The kinetic energy equation becomes 

1 



> — -^Ax 
h dt 



Ax 



4 



2 
Am 



( M; -+ M/+1 )Dj +| -D 



7+J 



Ax 



Ax 



Au i + l . 4 



Ax 



Ax 



The artificial dissipation terms in the flux lead to dissipation of kinetic energy since the 
first term on the right is negative while the remaining terms are consistent with the con- 
tinuous equation fl3.1) . Thus stability in the sense that the kinetic energy does not grow 
spuriously is maintained by the choice of the scalar dissipation flux. 



5.2 Entropy condition 

The dissipation in the energy flux will be determined by satisfying entropy condition 
approximately or exactly as before. The natural choice is to take D e = AE but this does 
not allow us to show entropy stability. Instead, the dissipation in the energy equation is 
taken to be of the form 

°' = { 2(^lj| + \"' U '» } A " +P + 2fFT) A < X/ « 
The entropy equation including the dissipation term is 

E^WEo( A * 3 )=-E 

; ; / 

\ 

— -Y"A, , i Av, , i D, , 1 
9 z r' /+2 1+2 l+i 

l 

The first term on the right is consistent with the continuous entropy equation while the 
second term is due to the dissipative flux. The contribution from the dissipation terms 



'Au 



'AT, 



Ax 



RTjT j+1 



i+l 
Ax 



Ax 
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can be calculated as follows. 



AvD = 



Ap 1 Aj6 



L P 7-1/3 
+2(pAu+uAfi) (pAu + uAp) 



fO(Ax 3 )-2 j 6uAi/-i/ 2 A J 6 



Ap 



-2A(j8) 
(Ap) 2 , 



2( 7 -l)/3 2 



\--UjUj + i > Ap+p uAu 



2pP(Auy + 



(7-1)^ 



2(7-1) 
(A/5) 2 +0(Ax 4 ) 



iyA(l/^) 



>o 

The leading order terms in the last equation are positive and hence lead to entropy gen- 
eration; the scheme satisfies the entropy condition in the limit of Ax — > 0. 

The scalar dissipation terms can also be constructed to be exactly entropy dissipative 
which can be used in combination with the entropy conservative flux. The dissipation 
in the mass and momentum fluxes are chosen as before which leads to the same kinetic 
energy stability as above, while the dissipation in the energy flux is taken to be 



D e = 



1 



2(7-1)0 2 



\--UjUj+i > Ap+p uAu + 



2(7-1) 



A(l//S) 



The only difference from the previous case is that the logarithmic average of jS has been 
used instead of the arithmetic average. Then the entropy production from the dissipative 
terms is 

(Ap) 2 



Av D = 



+2pP{AuY 



(7-l)/5;P>l 



(Apy >o 



and is non-negative which satisfies the entropy condition. We note that the dissipation is 
triggered if any one of the state variables p,u, /3 is non-uniform. The physical dissipation 
on the other hand acts only in the presence of velocity and /or temperature gradients, 
and is unaffected by density gradients. 



5.3 Summary of scalar dissipation flux 

The second order kinetic energy preserving and approximately entropy consistent scalar 
dissipation terms are given by 



D p 1 



DJ +1 _ = A(pu) j+l =u jH Ap j+l +p j+1 _Au j+1 _ 



D 
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and 

/ 'V 

(5.1) 




is the maximum wave speed at the interface. If we replace the arithmetic average jSy + i 
with the logarithmic average p^i in the above equations, then we obtain the exactly en- 
tropy consistent dissipation terms. Scalar artificial dissipation reduces the scheme to first 
order accuracy everywhere, so we can tune the dissipation to switch on only near shocks. 
Pressure gradients are a good indicator of the presence of shocks. However, lack of any 
dissipation in smooth regions can lead to instabilities. Hence, as in the JST scheme [13|, 
we will add a blend of second order and fourth order dissipation by replacing the differ- 
ence terms in the above equation, e.g., 

Ap /+ i ->e£\ (p j+ i-pj)-e®i (p j+2 -3pj+i+3pj-p h i) (5.2) 

where and e( 4 ) are adapted to the flow, with similar expressions for Aw and AT. Define 

\p h i-2pj+p j+1 \ 

v i = i 1 r\ I — r, v j+ 1 = max (v,-,V;+i 

1 \p h i+2pj+pj + i\ /+a 1 1 

e (2 \ =mm(l, K Wv i+l ), e {4 \ =max(0,K( 4 )-e (2) 1 ) 

]+ 2 I'l ]+2 1+2 

with k( 2 ) > and > 0. In smooth regions of flow, = 0(Ax 2 ), =0(1) and D = 
0(Ax 3 ), while near shocks =0(1), e^ 4 ' =0 and D = 0(Ax). Thus the scheme is second 
order accurate in smooth regions of the flow and if the mesh is well resolved, the scheme 
is second order accurate everywhere. For a very well resolved flow, the fourth order 
dissipation should be sufficient to stabilize the scheme as we already have consistency 
of kinetic energy and entropy condition which imparts some non-linear stability to the 
numerical scheme. 



6 Matrix dissipation flux 

Scalar dissipation schemes perform satisfactorily for weak shocks. For problems with 
strong shocks, one needs to add greater amount of scalar dissipation which unfortunately 
smears the contact discontinuities. An alternate way to add dissipation is via a matrix 
operator which leads to better shock capturing schemes. The classic example is the Roe 
scheme |20| which is based on a linearization of the non-linear conservation law about 
some average state. The numerical flux of the Roe scheme has the form 




19 



where R is the matrix of right eigenvectors and A is the diagonal matrix containing the 
eigenvalues, both of which are evaluated at the Roe average state. 



R- 



1 

u — a 
H—ua 



1 

u+a 
H+ua 



A| = |A| Koe = diag[|«-<i| / \u\, \u+a\] 



The dissipative flux can also be written in terms of the jumps in the entropy variables. By 
a linearization we can write Au = u v Av and due to a theorem of Barth [1 |, there exists a 
scaling of the eigenvectors R— >R such that u v =RR T . Then a Roe-type flux can be written 
by using jumps in the entropy variables instead of the conserved variables as 



i(f,+f /+1 )-i« /+5 |A ;+ ,|*7;,Av /+s 

If we want to use the eigenvectors in the original form R, then the flux can be written as 



h+\ 



-R f , i|A,, i|S,., iR , iAv,, i 

2 1+2 1 1+2 1 1+2 1+2 1+2 



2(f;+f/+0- 

where the matrix S provides the appropriate scaling and is given by 



(6.1) 



S = diag 



P_ (7~l)p 

2 T ' y 



_P_ 

'2 7 



In fact the matrix S is chosen so that R _1 du = SR T dv which also motivates the form of 
the flux given by equation ( |6.1| . Following this approach we can add a matrix dissipation 
flux to our kinetic energy and entropy consistent/ conservative centered flux f * to obtain 
the following dissipative numerical flux function 



f l+2 



f* 

1+1 



-R. , i |A - , i |S ; ,iR. iAv - , i 

2 l+2 l l+2 l 1+2 1+2 1+2 



(6.2) 



The above numerical flux together with | A| = | A| Roe will be called the KEP-ES flux; when 
the approximately entropy consistent flux is used for the central part, it will be explicitly 
qualified as (AC). We note that the dissipation matrix Q = R|A|SR T is positive definite 
which allows us to derive the entropy inequality as follows. By taking the dot product of 
the semi-discrete finite volume scheme for Euler equations with vy we obtain the entropy 
equation 

dUi 1 



Ax- 



df , Vr{ i+i 2 
Define the numerical entropy flux to be 
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which is a consistent flux, then the cell entropy equation can be written as 

dUj _ _ _,. „ 1 

~4 



Ax- 



-F. i- 

J 2 



-0(Ax 3 



-Av 



iAv, 

2 / 



<0 



dt ' ~J- 

which shows that the dissipation terms in the flux lead to entropy generation. The 
O(Ax) 3 terms are not present if we use the exactly entropy conservative flux for the cen- 
tered flux f * and we then obtain strict entropy stability for any mesh size Ax. 



6.1 Resolution of stationary contact discontinuity 



Contact discontinuities occur because of the presence of linearly degenerate eigenvector 
fields. Because of their linear nature, the accurate resolution of contact waves is more 
difficult than shocks which have an inherent steepening mechanism due to their non- 
linear nature. The original Roe scheme exactly preserves stationary contact waves. This 
is a desirable property since it affects the accuracy with which shear layers and boundary 
layers are computed. The new flux functions with entropy variable based dissipation 
will be able to resolve stationary contacts exactly if we evaluate the enthalpy H occuring 
in the eigenvector matrix R in an appropriate manner. Consider an initial condition such 
that 

( PI k <) 

fr k>j' 

This initial discontinuity satisfies the RH jump conditions with zero speed and hence it 
represents a stationary solution. The numerical scheme will preserve this solution exactly 
if 

"ol 

. = = ij = f j+ i = f . + 3 = f ; - +1 : 



Pk = 



u k = 0, pk = P 



f ,+r 



For the above initial condition, the centered flux f * already satisfies the above condition, 
i.e., f*^i = [0, p, 0] T . Hence we need to have Qy + i Avy + i =0 in order to preserve the contact 



;+ 



wave. The dissipation matrix Q requires some average values of u, a, H and p which 
will be denoted by u- + i, etc. For the above initial condition, any consistent averaging 



would yield u- + 1 =0. We will compute the enthalpy as H- + 1 = a 2 . + 1 / (7 — 1 ) + u 2 . + j / 2. The 



dissipation matrix then takes the form 



so that 
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This is zero provided the term inside the brackets is zero. We will show that this is the 
case if we choose 



7 



2 ^ 



(6.3) 



Then 



7" 



Ap. 



1 A/3 ; , 
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+1 
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where the last equality follows due to the chosen initial condition. Thus exact resolution 
of stationary contacts is possible provided we evaluate the sound speed a- + \ using the 
logarithmic average of f> as in equation (6.3 1. The other quantities Py + i, u.- + i can be eval- 
uated using arithmetic or logarithmic averaging. Since the velocity already appears in 
the central flux f * as an arithmetic average we can use the same choice in the dissipation 
matrix also. The density appears in the approximately consistent flux as an arithmetic 
average and in the exactly conservative flux as the logarithmic average; we will use the 
corresponding approximations in the dissipation matrix which avoids additional com- 
putations. 



6.2 Kinetic energy stability 

Consider the finite volume scheme with entropy consistent/ conservative flux together 
with entropy variables based matrix dissipation. The semi-discrete kinetic energy equa- 
tion for this scheme is 



£^r Ax = -2E(";- M ; + i)(2( M / +M /+i)(Q Av ); + r (QAv) rH 



Au m _ 4 



u-.i, -1, 
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The second term is consistent with the exact kinetic energy equation while the first term 
is due to the additional dissipation in the flux. This term will be dissipative for the kinetic 
energy provided 



u j+ i, -1, 



i >0 

J1~2 



(6.4) 
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Now dropping subscripts and overbars, we have 

ap 
27 



[u, -1, 0]QAv = ^[(|Ai|-|A 3 |)At> 1 -(a(|Ai| + |A3|)+u(-|A 1 | + |A 3 |))At>2 



-(a«(|Ai| + |A 3 |)+H(-|Ai| + |A3|))A©3] 

Note that since Av\ contains Ap whereas AV2, Av$ do not, we can satisfy condition ( |6.4| 
only if I Ai I = I A3 1 = A for some A > 0. Then the kinetic energy equation becomes 



dK 1 
T —r-&x = - — Ya 1 l ip., 16., iA,, 1 (Aw,., i) 2 + Y" 

7 ' 7 7 



A M/+ i_ 4 /Am + i\ 2 



Ax 



The first term on the right is dissipative and hence the scheme is stable for the kinetic 
energy. The Roe-type dissipation does not satisfy the required condition since \\ = u — 
a and A 3 = u+a. The condition that the first and third eigenvalues should be equal is 
satisfied if we choose the Rusanov form of the dissipation which corresponds to 

\A\ = \A\ Rus = Xl, X=\u\+a (6.5) 

However this scheme is very dissipative and leads to excessive smearing of shocks, con- 
tact discontinuities and shear layers. Since only the first and third eigenvalues have to be 
equal to have dissipation of kinetic energy, another obvious choice is 

|A| = |A| K£S = diag[ A, \u\, A ], A=max(|w — a|,|w+a|) = |tt|+a (6.6) 

which resolves stationary contacts exactly. However steady shocks are smeared over 
many cells due to excessive dissipation. The robustness of the scheme due to its kinetic 
energy and entropy stability makes it attractive when used in combination with high 
resolution schemes like WENO or discontinuous Galerkin methods. The new entropy 
conservative flux toegther with the above forms of dissipation will be refered to as KEP- 
ES(Rus) and KEP-ES(KES) respectively. 



7 Monotone resolution of shocks 

An important desirable property of a numerical scheme for hyperbolic problems is that 
it should yield monotone discontinuous solutions. Central schemes yield highly oscilla- 
tory solutions near shocks and other discontinuities which are physically incorrect and 
can even lead to negative density and pressure. Upwind schemes which are based on 
wave propagation ideas and central schemes with explicitly added dissipation are able 
to yield non-oscillatory solutions. For scalar conservation laws, the TVD condition pro- 
vides a complete solution for the design of accurate, non-oscillatory schemes. However 
the TVD property does not carry over to non-linear systems like the Euler equations and 
there is a lack of theoretical basis for the design of non-oscillatory schemes for systems 
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of conservation laws. It is hoped that the entropy condition can be a useful criterion for 
the design of non-oscillatory schemes. However this condition only specifies that there 
must be non-zero dissipation which leads to some entropy production but does not say 
how much entropy generation is necessary to yield monotone solutions. The matrix dis- 
sipation flux discussed in the previous sections may not yield monotone shock solutions 
even though they are entropy stable. In order to study this numerically we follow Ismail 
and Roe [ 10 1 and consider a problem with a stationary shock solution. This corresponds 
to a Riemann problem with the following left and right states 



Pl = l, Ui = l, 



Pi 



and 



7-1 



(7+l)Mf 7+1 
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2yMf 7-1 
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The Mach number on the left M/ can be changed to produce shocks of different strengths 
and following |l0) we consider M/ = 1.5, 4 and 20. All computations are performed using 
a CFL=0.1 and a three stage Runge-Kutta scheme. The density obtained from the different 
schemes is shown on the left of figure jl]). The ROE-ES flux and the new approximately 
entropy consistent flux gives pre-shock oscillations in the density at all Mach numbers 
while the new entropy conservative flux yields oscillations at the lower Mach number 
but is monotone at higher Mach numbers. Ismail and Roe attribute the non-monotonicity 
to insufficient entropy production and hence insufficient dissipation at the shock. For a 
weak shock, the entropy production is O(Ap) 3 whereas the entropy production due to 
the dissipative flux is O(Ap) 2 . Hence they propose modifying the acoustic eigenvalues 
so that the eigenvalue matrix |A| becomes 

|A| = |A| EC1 =diag[|«-a|+^|AAi| / \u\, \u+a\+^\^\\ (7.1) 

where AAi and AA3 are the jumps in the corresponding eigenvalues across the cell face. 
If j6 = g then the entropy production across a weak shock due to matrix dissipation corre- 
sponds to the correct entropy production of O(Ap) 3 as shown in [ 10 1 . With this modifica- 
tion of the eigenvalues, the entropy conservative Roe flux is called ROE-ECl flux where 
EC stands for "entropy consistent". We will adopt the same eigenvalue modifications 



as in equation (7.1 1 in combination with the new entropy conservative fluxes refered to 
as KEP-EC1 and repeat the computation of the stationary shock problem. The solutions 
shown on the right of figure jl} indicate that the ROE-ECl flux is able to give mono- 
tone solution at the lower Mach number but still produces some oscillations at the higher 
Mach numbers. The new approximately entropy consistent flux produces oscillations at 
all Mach number^] while the new entropy conservative flux yields monotone solutions 



* These oscillations disappear if we use /3 = 1. 
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at all Mach numbers. This shows that exact entropy conservation does yield better solu- 
tions and is to be prefered. In order to improve the monotonicity of the ROE-EC 1 flux, the 
acoustic eigenvalues can be further augmented leading to the ROE-EC2 flux as discussed 
in 1 10 1. However with the new entropy conservative flux, we find that there is no need 
to increase the eigenvalues for higher Mach numbers and the basic Roe-type dissipation 
yields monotone solutions. The behaviour of the new flux is consistent with the entropy 
production analysis of p0| which is valid only for weak shocks. 

We also compute the problem using the entropy dissipation that leads to kinetic en- 
ergy stable scheme. Note that we do not use augmented eigenvalues as in the ROE-EC1 
scheme but the first and third eigenvalues are taken to be equal to the maximum eigen- 
value while the second eigenvalue is unchanged as in equation ( |6.6| . The solution for 
Mach = 20 is shown in figure ([2]); the first order scheme is highly diffused due to the in- 
creased dissipation and the scheme is non-oscillatory. If we use the linear reconstruction 
scheme with minmod limiter, then a sharper shock resolution is obtained which is again 
non-oscillatory. Such kinetic energy and entropy stable schemes would be attractive for 
use with high order accurate schemes like WENO and discontinuous Galerkin methods. 



8 Shock instability, carbuncle and a hybrid scheme 

Numerical schemes based on exact or approximate Riemann solvers suffer from a variety 



of pathological behaviours 1 19 1. The entropy violating shock can be overcome by ensur- 
ing that the numerical scheme is consistent with the entropy condition. Other problems 
like the shock instability and carbuncle do not have a proper rational fix available. 

Consider the stationary shock problem as discussed in the previous section. If we set 
the initial conditions to exactly correspond to the discontinuity, then all the schemes yield 
stable, stationary solutions for which the residual converges to machine zero. However 
the Roe scheme and many other schemes are known to be unstable if the initial condition 
contains an intermediate point at the shock location. In order to induce the instability, 
the mass flux at the right boundary, which is an outflow boundary, is fixed to be one, 
while the momentum and energy fluxes are computed by zero gradient condition; this 
amounts to not updating the momentum and energy in the last cell. Small disturbances 
produced at the shock due to the intermediate point are propagated downstream, get 
reflected from the outflow boundary back into the shock and further disturb the shock 
location leading to a limit cycle oscillation. The new KEP-ES schemes do not suffer from 
the shock instability problem though the solutions are not monotone. The addition of 
extra dissipation in the KEP-EC1 scheme leads to the appearance of instability at higher 
Mach numbers. 

It is thought that numerical schemes which suffer from this instability will also be 
prone to the carbuncle phenomenon which is found in flows with strong shocks over 
blunt bodies. However there is no guarantee that a scheme which is stable for the 1-d 
shock problem will not produce the carbuncle phenomenon. The ROE-EC2 flux which 
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Figure 1: Stationary shock problem: N = 24 cells 
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Figure 2: Stationary shock problem: N = 24 cells 



has increased dissipation avoids the 1-D shock instability but can still suffer from the 



carbuncle problem in multi-dimensions |11|. We show in the later sections that the KEP- 



EC1 scheme also produces carbuncle effect on certain types of meshes where the shock is 
aligned with the mesh. Increasing the value of /3 in equation (7.1 1 even upto /3 = 1 does 
not seem to eliminate this problem. The kinetic energy stable scheme given by equa- 
tion ( |6.6| which has more dissipation in the acoustic waves also produces the carbuncle 
but the Rusanov scheme of equation ( |6.5| does not produce the carbuncle. In fact what 
we observe is that all the schemes which resolve grid aligned stationary contacts exactly 
seem to suffer from the carbuncle effect which is consistent with what is noticed in other 
schemes in the literature. The usual fix in such cases is to increase the amount of dissipa- 
tion in the numerical scheme which however causes poor resolution of boundary layers. 
There is also the idea of switching the numerical scheme to a more dissipative one only 
near shocks and using a high resolution Riemann solver type scheme in smooth parts of 
the flow [19|. In the framework of the entropy conservative /stable scheme as discussed 
in this paper, we have the freedom in designing the eigenvalues which essentially control 
the amount of dissipation in the scheme. Since the Rusanov scheme is free of carbuncles, 
we propose a blending of the usual Roe scheme with the Rusanov scheme, i.e., the matrix 
| A | appearing in the entropy dissipation flux is of the form 



|A| = |A|% fc = (1-0) | A\ Roe +(p\A\ 



Rus 



where the switching function cp is based on the pressure jump 

i 

Ap 



2p 



Note that 0<(p<l; for a strong shock i^l and the scheme is close to the Rusanov scheme, 
while near a contact discontinuity (p ~ and the scheme is close to the more accurate Roe 
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scheme. This hybrid scheme does not produce the 1-D shock instability problem. In 
the results section we show that this scheme avoids the carbuncle effect and still gives 
good resolution of boundary layers and shear layers. The new entropy conservative flux 
together with the above hybrid dissipation will be refered to as KEP-ES(Hyb) scheme. 
We also remark that these modifications of the dissipation term still retain the entropy 
stability property of the scheme. 



9 2-D NS equations and finite volume method 



The two dimensional Navier-Stokes equations in conservation form can be written as 
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with u = {u\,U2) being the Cartesian components of the fluid velocity and E = p/ (7 — 1)- 
p\ u\ 2 /2 is the total energy per unit volume, while the viscous fluxes are given by 
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where Ty is the shear stress tensor and c\i is the heat flux vector, for which we assume 
the Newtonian and Fourier constitutive laws, respectively. We will approximate these 
equations using a finite volume method on unstructured grids. Specifically we will use a 
vertex-based finite volume scheme in which a primal grid of triangles is used; the finite 
volumes are constructed around each vertex using either a median dual cell or a voronoi cell. 
In the median dual cell, the finite volumes are constructed by joining the cell centroids 
to the edge mid-points while in the voronoi dual, the circumcenter of the triangles are 
joined to form the finite volumes. If the triangle is obtuse angled so that the circumcenter 
lies outside the triangle, then the mid-point of the largest side is used instead of the 
circumcenter. 

Let Aj be the area of the f'th vertex/ finite volume and let N(i) denote the set of neig- 
bouring finite volumes which share a common boundary with A/. We also use the set 
C(z') which is the set of all triangles T containing the z'th vertex. The semi-discrete finite 
volume scheme is given by 

U ' jGN(i) TGC(i) 
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where f« is the flux from the z'th cell into the j'th cell across their common boundary 
whose length is As,,-. The discretization of viscous fluxes uses a Pi finite element approach 
on triangular grids which leads to only a nearest neighour stencil (edge connected neigh- 
bours) and hence is very compact. The flux gf is the contribution of the dissipative flux 
to vertex i coming from the portion of the boundary of A\ lying inside triangle T and As ; T 
is the length of this boundary. For a boundary finite volume, there will be flux contri- 
butions from the boundary edges also. We will approximate this flux using the entropy 
conservative flux together with some dissipation which takes the form 

The entropy conservative flux f* in the direction of the unit normal vector n = (wi,^) to 
the cell face is given by 

p(u-n) 
prii+Uif*'P 
pn 2 + U2f*' p 
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where p, ft are the logarithmic averages and p is given by equation (4.11 
eigenvectors R is given by 



The matrix of 
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and the scaling matrix S is 
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The eigenvalue matrix for the basic Roe scheme (ES) is given by 

|A| = |A| Roe = diag [|w-n — a\ \u-n\ \u-n\ \u-n+a\] 

and the modifications for the EC1, KES, Rus and Hyb schemes are obvious from their 1-D 
counterparts. Finally the entropy variables are given by 



7-1 



-p\u\ 2 2/3wi 2/3m 2 -2j6 



Second order accuracy is achieved by a MUSCL type reconstruction using nodal gra- 
dients of primitive variables p,ui,ii2,T which are calculated using Green-Gauss theo- 
rem while time-stepping is done using a 3-stage strong stability preserving Runge-Kutta 
scheme [23 [. Triangular grids used in the computations are generated with GMSH [6] 
and Delaundo 1 17 [. 
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Figure 3: Modified Sod problem: Density at time f = 0.2. The figures on the right show a zoomed view of the 
rarefaction region. 



10 Numerical results 

In all the test problems, we consider an ideal gas with 7 = 1.4 except when indicated 
otherwise. 

10.1 Modified Sod test case 

This is a shock tube problem with the left state being (p,u,p) = (1.0,0.75,1-0) and the right 
state being (p,u,p) = (0.125,0.0,0.1). All computations are made with N = 100 eels and a 
CFL=0.4 upto a final time of t = 0.2 units. The original Roe scheme gives entropy violat- 
ing jump in the expansion region where the flow becomes sonic as shown in figure (pV 
The reason for this is the vanishing of one eigenvalue at the sonic point which leads to 
insufficient dissipation. All the entropy consistent schemes including the new approxi- 
mately entropy consistent scheme together with entropy variable based dissipation give 
solutions which do not suffer from this problem as shown in figure (p). The zoomed view 
of the figure on the right shows that in the sonic region of the expansion fan all of them 
behave in a very similar manner. Note that even the entropy stable schemes have a van- 
ishing eigenvalue in the expansion fan. However due to the entropy conservative nature 
of the central part of the flux f *, they do not give rise to the entropy violating shock unlike 
the original Roe scheme which is not entropy conservative. 
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Figure 4: Stationary contact 

10.2 Stationary contact 

We consider inviscid Riemann problem with left state (pi,ui,pi) = (10,0,1) and right state 
(p r ,u r ,p r ) = (1,0,1) whose solution is a stationary contact wave. The problem is solved 
with N = 26 cells and the solution is shown in figure (Q. The condition for the exact 
resolution of stationary contacts is the correct computation of the enthalpy in the dissi- 
pation matrix. For the approximately entropy consistent flux (AC), we use the arithmetic 
average of f> to compute the enthalpy which does not satisfy the exact contact resolution 
property. For the exactly conservative flux, we use the logarithmic average of f5 which 
should resolve the contact wave exactly. The exactly conservative flux is also used in 
combination with the kinetic energy stable scheme KEP-ES(KES) which should resolve 
the contact discontinuity exactly. The results of these three schemes are consistent with 
the theoretical prediction as seen in the figure These results show that contact waves 
can be highly dissipated if the numerical flux is not able to resolve contacts exactly. 

10.3 Sod test case 

Inviscid case: This is also a shock tube problem with the left state being (p,u,p) = (1.0,0.0,1.0) 
and the right state being (p,u,p) = (0. 125,0.0,0.1). All computations are made with N=100 
cells and a CFL=0.4 and a 3-stage Runge-Kutta scheme upto a final time of f = 0.2 units. 
We compute the solution using entropy consistent scheme with entropy dissipation and 
using MUSCL scheme and minmod limiter. The solution shown in figure |5j indicates 
that the sharp resolution of contact and shocks can be achieved with the new schemes. 
The shock is resolved within two cells and the contact is resolved with about four cells. 
It can be seen that all the schemes give essentially the same solution on this problem. 
Viscous case: The same problem is solved using Navier-Stokes equations at a Reynolds 
number of 2000 based on the sound speed of the left state and the length of the domain 
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using 500 cells and a CFL number of 0.1. We first use the central fluxes without any dissi- 
pation; the solutions are obtained with the kinetic energy preserving scheme of Jameson 
(KEP), Roe's entropy conservative scheme ROE-ES and the new KEP-ES scheme. The 
entropy s = ln(pp~ 7 ) is shown in figure |6^) in a zoomed section in the expansion region. 
The density shown in figure |6j?) also oscillates in the expansion region but the oscilla- 
tions are too small to be seen on the scale of the figure. It is found that all the central 
schemes produce some oscillations in density and pressure which is also reflected in the 
entropy while there are no oscillations in velocity and temperature. Reducing the CFL 
number does not eliminate these oscillations. These oscillations originate in the discon- 
tinuous initial condition, particularly density and pressure, and propagate upstream. The 
KEP scheme is seen to produce smaller oscillations compared to the entropy conservative 
schemes. This indicates that the KEP scheme has some inherent dissipation which is able 
to damp the oscillations in density but it is not able to eliminate them completely. The 
entropy conservative schemes do not have this additional dissipation probably because 
of their entropy conservative nature. Note that the physical viscosity which acts due to 
the velocity and temperature gradients is ineffective in damping oscillations in density. 
This problem is present even with implicit kinetic energy preserving schemes and in |24) 
the authors modify the pressure flux to be biased instead of centered; this adds some dis- 
sipation which becomes active if there are oscillations in pressure. In the present work, 
we also perform computations using the KEP-ES scheme but adding only fourth order 
scalar dissipation (D4) with k' 4 ' = j^j. As shown in figure roL the fourth order dissipation 
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(a) (b) 

Figure 6: Sod problem using Navier-Stokes equations: Solution at time t = 0.2 using 500 cells (a) Zoomed view 
of entropy (b) Density 



is able to damp the density oscillations without affecting the accuracy of the solution in 
other regions. In particular the zoomed view in the shock region shows that the fourth 
order dissipation does not spoil the accuracy of the solution even inside the shock region. 



10.4 NS shock structure 

We compute the shock structure using the Navier-Stokes equations with the help of the 
kinetic energy preserving and entropy conservative flux together with scalar artificial 
dissipation. In the artifificial dissipation we use = | and = ^. The parameters 
defining the problem are: Mach number ahead of the shock is M\ = 1.5, 7 = 5/3, Pr = 2/3 
while the viscosity law is given by \i = ji\{T /T\)°- % where the subscript "1" denotes pre- 
shock conditions and \i\ =0.0005. Figure 0-((9]l shows the solutions on N = 50,100,200 
cells. On the coarse mesh, the stress and heat flux cannot be computed accurately since 
there are too few points inside the shock region, but the solution is still non-oscillatory. 
On the finer meshes, the scheme is able to compute the shock structure with good accu- 
racy. Using N=200 cells, we compute the solution with only the fourth order dissipation, 



K ^ = 200 anc ^ * ne solution is shown in figure ( 10 >. We are thus able to obtain accurate so- 
lutions on fine meshes for the Navier-Stokes equations using the central kinetic energy 
and entropy conservative scheme. The fourth order dissipation helps to damp the oscil- 
lations in density/ pressure which are created due to the initial discontinuity as discussed 
in the Sod test case. 




re 8: IMS shock structure: N = 100 cells, KEP-ES(SD) flux with second and fourth order dissipation 




Figure 10: NS shock structure: N = 200 cells, KEP-ES(SD) flux with fourth order dissipation 
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(a) (b) (c) 

Figure 11: Grid used for supersonic cylinder problem: (a) Primal grid, (b) median dual grid and (c) voronoi 
dual grid 



10.5 Supersonic flow past cylinder 

Consider the inviscid supersonic flow over a semi-cylinder; the primal triangular grid 
and corresponding median and voronoi dual meshes are shown in figure ( |TTj |. The 
voronoi cells lead to nearly structured type grids and can thus lead to carbuncle prob- 
lem since the shock will be aligned with the cell faces in a more exact way than for the 
median dual cells. For a free-stream Mach number of 2 on the voronoi cells, the solution 
of the KEP-EC1 scheme is shown in figure ( [l2] | and there is no problem of carbuncle ef- 
fect. For a Mach number of 20 and using the median dual grid the KEP-EC1 scheme is 
able to give carbuncle free solutions with good shock resolution as shown in figure p3| . 
However the same scheme gives rise to carbuncle effect on the voronoi mesh as seen in 



figure ( 14 1). The carbuncle is also seen in the kinetic energy stable scheme KEP-ES(KEPS) 



in figure (14 ?). These two schemes resolve stationary contacts exactly and they also ex- 



hibit the carbuncle effect. The Rusanov and the new hybrid schemes do not give the 



carbuncle problem as seen in figures ( 14 :) and (|14p) respectively. By looking at the shock 
thickness, it can also be observed that the hybrid scheme is less dissipative than the Ru- 
sanov scheme, and this is further confirmed in the boundary layer test case. We also 
remark that in our numerical experiments the KEP-ES scheme (i.e. without any augmen- 
tation of the eigenvalues) does not produce the carbuncle effect even on the voronoi cells 
but the solution is not monotone as already seen in the 1-D test cases. 

10.6 Transonic flow past NACA-0012 airfoil 

This is a standard test case for inviscid aerodynamic problems and involves a symmetric 
NACA-0012 airfoil placed in a freestream Mach number of 0.85 and angle of attack (AOA) 
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(a) (b) (c) (d) 

Figure 14: Supersonic cylinder, Mach=20: Voronoi dual grid, density contours (a) KEP-EC1 (b) KEP-ES(KES) 
(c) KEP-ES(RUS) (d) KEP-ES(Hyb) 
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Inlet portion Flat plate 

Figure 16: Primal grid for laminar flat plate boundary layer problem 



of 2 degrees. The flow develops shocks both on the upper and lower airfoil surfaces. 
We compute this flow on a triangular grid containing 180 points on the airfoil surface 
and 20 points on the farfield boundary which is a circle, with a total of 6402 vertices. 
Higher order accuracy is achieved using MUSCL type reconstruction and van Albada 
limiter. The solution using the original ROE scheme, the new KEP-EC1 and KEP-ES(Hyb) 
schemes are shown in figures p5| . The results with ROE and KEP-EC1 schemes are 
almost identical, while the KEP-ES(Hyb) scheme which has increased dissipation shows 
less sharp resolution of shocks. For such transonic flows, it is preferable to use the KEP- 
EC1 scheme which is accurate and satisfied entropy condition. 



10.7 Laminar flat-plate boundary layer 

This problem corresponds to viscous flow over a flat plate which leads to the develop- 
ment of a boundary layer near the plate surface. The Reynolds number corresponding to 
the plate length is 10 5 while the Mach number of the incoming flow is taken to be 0.1. The 
computation domain is rectangular as shown in figure ( 16 1 which also shows the primal 
triangular grid used for the computations. There is an initial inlet portion of the domain 
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Figure 17: Laminar flat plate boundary layer: (a) Streamwise velocity, (b) Vertical velocity 



on which slip boundary condition is imposed followed by the no-slip boundary corre- 
sponding to the flat plate. Adiabatic conditions are used on the flat plate boundary. At 
the top and outlet, the free-stream pressure is specified while at the inlet the free-stream 
values are used together with the numerical flux function to compute the flux. We com- 
pare the numerical solution of the velocities with the Blasius semi-analytical solution in 
figure py} in the standard non-dimensional units. These results are taken on the vertical 
line through the center point of the plate. We have shown the numerical results using the 
Rusanov and hybrid form of the dissipation terms. The results from ROE and KEP-EC1 
schemes are almost similar to the KEP-ES(Hyb) scheme and are not shown here. It is clear 
that the new hybrid scheme is able to give accurate resolution of the boundary layer pro- 
file while the Rusanov scheme introduces too much dissipation. The difference is more 
dramatic in the vertical velocity profile which is more sensitive to numerical dissipation 
since the vertical velocity component is much weaker than the streamwise component. 



10.8 Step in wind tunnel 

This test case is described in [28 1 and involves inviscid supersonic flow past a step in a 
wind tunnel which is impulsively started. The initial Mach number is 3. The flow devel- 
ops several shocks which undergo reflections. A shock triple point intersection leads to 
the formation of a slip line. The corner on the step is a singular point and many numer- 
ical schemes develop spurious entropy at the corner which then leads to the formation 
of a Mach stem in the downstream direction. This problem is solved using the KEP-EC1 
and KEP-ES(Hyb) flux functions together with linear reconstruction of primitive vari- 
ables and minmax limiter p). The grid is adapted to be finer near the corner where the 
spacing is of 0(0.002) while the maximum spacing is of 0(0.01), the total number of grid 
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Figure 18: Forward step in wind tunnel at M = 3: Mach number contours, 50 equally spaced contours between 
and 4.8, (a) KEP-EC1 flux, (b) KEP-ES(Hyb) flux 



points being 64246. The Mach number contours at time t = 4 are shown in figure (18} 
using the KEP-EC1 and KEP-ES(Hyb) schemes and both of them are able to resolve the 
main features of the flow. The slip line is resolved equally well by both schemes and the 
KEP-ECl(Hyb) does not add any extra dissipation at the contact waves. The reflection of 
the shock from the botton wall is seen to be free of the spurious Mach stem. 

11 Summary and conclusions 

We have derived a new entropy conservative flux for the Euler equations which also 
preserves the kinetic energy. Due to their better consistency properties, such numerical 
fluxes could be attractive for direct numerical simulation of Navier-Stokes equations on 
highly resolved meshes. These fluxes are central in nature and hence lead to second order 
accurate schemes. But due to their central character, for coarse meshes or for problems 
in which the shock is not well resolved by the mesh, some dissipation has to be added to 
stabilize the scheme. We have constructed scalar and matrix dissipation schemes which 
preserve the stability properties. In particular, dissipation of kinetic energy is shown to 
occur if the eigenvalues in the dissipation matrix are chosen appropriately. The matrix 
dissipation schemes are based on Roe-type dissipation using entropy variables; how- 
ever even the first order schemes can yield oscillatory solutions even for weak shocks. 
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The eigenvalues can be augmented as proposed by Roe which leads to monotone resolu- 
tion of shocks even at hypersonic Mach numbers. Due to the entropy consistency of the 
central part of the numerical flux, all of the entropy consistent schemes avoid unphysical 
shocks at sonic points. However, like other schemes based on Riemann solvers and which 
resolve stationary contacts exactly the new scheme can also suffer from 1-D shock insta- 
bility and the multi-dimensional carbuncle problem at higher Mach numbers. A hybrid 
dissipation scheme which blends the Roe type dissipation with the Rusanov dissipation 
is proposed and shown to be free of carbuncle effect for the blunt body problem. The 
blending is achieved at the level of the eigenvalues or wave speeds which appear in the 
dissipation matrix. The hybrid scheme is also able to give good resolution of boundary 
layer flows which is an important requirement in the computation of viscous flows and 
is more accurate than the Rusanov type dissipation scheme. Several standard test cases 
that are commonly used for validation of compressible flows have been computed with 
the new class of schemes and shown to yield accurate solutions. 

Kinetic energy and entropy conservative schemes are attractive for DNS of shock- 
free compressible turbulent flows. It was shown in the Sod shock tube problem solved 
with Navier-Stokes equations without any artificial viscosity (Section 10.3} that the KEP 
scheme gives smaller entropy oscillations than the entropy conservative schemes, indi- 
cating that there is some inherent entropy dissipation in the KEP scheme which is not 
present in the entropy conservative schemes. However the KEP scheme does not com- 
pletely eliminate these oscillations and there is a need to add some dissipation to control 
the oscillations in density and pressure. It is preferable to use the kinetic energy and 
entropy conservative scheme since it does not add any implicit dissipation; explicit dissi- 
pation to control density oscillations can be added whose magnitude is precisely known 
and hence can be controlled. 
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